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Abstract 

Using simple methods of asymptotic analysis it is shown that for a 
quantum harmonic oscillator in n-th energy eigenstate the probability 
of tunneling into the classically forbidden region obeys an unexpected 
but simple asymptotic formula: the leading term is inversely propor¬ 
tional to the cube root of n. 
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1 Introduction 

Consider the quantum harmonic oscillator in dimensionless variables, where 
the Hamiltonian, in terms of position and momentum operators x,p, satis¬ 
fying [x,p] = h is given by the formula 

Its normalized eigenstates, when represented as square integrable functions 
of the position variable x are 

Hn are Hermite’s polynomials: Hn{x) = (—. The correspond¬ 
ing probability densities Pn{x) are then given by |V’n(x)p, i.e. 
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where 


( 1 ) 


1 


— 


7ri/22’^n!' 

so that the normalization condition holds: 


/' 


Pn{x)dx = 1 . 


The functions Pn{x) are symmetric with respect to the origin x = 0. The 
classical turning points are at x = Pyj2n + 1, thus the probability of quan¬ 
tum tunneling into the classically forbidden region are given by the formula 




n.tun 


= 2 


poo 

Jj2n+l 


dx — ‘2(XyiQ'f. 
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( 3 ) 


where 

/ \ 2 

Qn= dx. 

Jv/2H+T ^ ^ 

The probability distributions Pn{x) are usually represented (see, for exam¬ 
ple, [U p. 168]) as in Fig. [T] 

The tunneling probabilities Pn,tun are rarely discussed in textbooks on 
quantum mechanics. Sometimes they are discussed in quantum chemistry 
textbooks. For instance in [3 p. 92] the first three are calculated (though 
the second one with an error: it should be 0.1116 instead of 0.116). In [8l p. 
66-67; p. 98-99 in 2nd ed (2014)] Po,tun is calculated with a comment: “The 
probability of being found in classically forbidden regions decreases quickly 
with increasing n, and vanishes entirely as n approaches infinity.” Yet the 
question of “how quickly” is left open there. To answer this question we will 
derive the asymptotic behavior, first by using a heuristic method based on 
a simple formula by Szego, then, in a rigorous way, using another formula 
due to Olver. Both methods lead to the same final result. The problem 
that we solve is closely related to a more general problem of large quantum 
numbers and the correspondence principle. As can be seen from Refs. la^s], 
even for such a simple quantum system as the harmonic oscillator, there are 
surprises (a more general discussion can be found in the monograph m)- 
Our result showing a very slow approach to zero of Pn,tun for large n may 
be considered as another surprise. 

In section [2] we derive the leading term formula using heuristic reasoning. 
In section [3] we use mathematically rigorous arguments based on dominated 
convergence theorem, while controlling the order of terms. In section l4] we 
derive the second order correction to the leading term. The proof of the 
main lemma used in this paper is given in the Appendix. Our theoretical 
results are also verified numerically. 
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n = 40 



Figure 1: Probability distribution P^oix), with the classical turning points 
scaled to x = ±1, and the curve Pn,tun for n from 0 to 40. The black [/-shaped 
curve represents classical probability distribution Pciass{x) = l/{'Ky/l — ) 

for an oscillator with known energy and unknown phase. Shadowed tails 
beyond the classical turning points are responsible for non-zero quantum 
tunneling probabilities. 


2 Derivation of the asymptotic tunneling proba¬ 
bility formula 


We are interested in the approximate behavior of Pn,tun, as a function of n, 
for “large” n. To this end we first use the asymptotic formula by Szego P P- 
201, (8.22.14)] ( valid in the transition regions ) 


-^HJx) = 


1 n 

3322 ’ 


■3(n!)i 


7r4ni2 


A(t} 1 1 -|- 0(n 3 


(4) 


x = ( 2 n + 1)^/2 _ 2-1723-1/3^-1/6^^ 

The function A(x) used by Szego relates to the Airy functions Ai(x) by the 
following formula (cf. [TOl p.l94]): 

A(t) = 


3 





































































( 5 ) 


Thus Eq. 0] can be written 


where 


s = - V2n + 1). 


The above asymptotic formula is proven to be valid for t bounded. In the 
calculation below I will use this result for all s > 0, because the main 
contributions to integral come from a very small neighborhood of s = 0, in 
particular when n is large in (jSj). 

Since ^ we can calculate now Qn dehned in Eq. [3] as 



We get 



The integral of the squared Airy function is known (cf. e.g. |12l Eq. (3.75), 
p. 50]): 



Therefore (cf. Eqs. (I2|),([T])) 


(7) 



1 


32/3 r 2 [i/ 3 ] 


_ 7 _ ri-^/3 

32/3r2[l/3] 


( 8 ) 


Thus, numerically. 


0.133975 


(9) 



As it can be seen from Eig. [21 the asymptotic formula ([9]) seems to agree 
with the numerical calculations for n up to 612. Yet in its derivation we have 
assumed the validity of Szego’s asymptotic formula in the unbounded inter¬ 
val. Can such an extension be mathematically justified? We will address 
this question in the next section. 

^Another use of the Airy function in a discussion of the harmonic oscillator wave func¬ 
tions can be found, for instance, in HU Fig. 5 - 23 ], in the section on WKB approximation. 
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Figure 2: Comparison of the tunneling probabilities calculated from the 
exact formula with the curve given by Eq. Q, starting from n = 5. 


3 Using Giver’s formula 

Olver m p. 403] gives the following formula valid for x > 1: 

e~^Hnii'x) = Cn ^ ^2*^ ^ (Ai(z^5C) + , (10) 

where 

. - 1 v'^ (3v^-l) 

v = \/2n +1, c,i = (27r)2e v s ^ (11) 

/ o o \ 2/3 

C = \^-^x{x^ - 1)^'^^ - ^ cosh"^(x) J , (12) 

and 

|e(x)| < 1.36 (exp(0.09z/"^) - l) Ai(z 2 ^/\). (13) 

Remark 1. ; The function cosh.~^ {x) in the above formula denotes the posi¬ 
tive branch of the inverse function o/cosh. For general aspects of asymptotic 
methods for integrals see m where, in particular, Section 23.4 provides 
a complete description of the Airy-type expansion of Hermite polynomials. 
Notice also that = l/(2n + 1) = (l/(2n))(l — l/(2n)) + 0{n~^). 
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We will need Eq. (fTO|) squared. Using Eqs. (fTO]) - (fT^ we write it in the 
form 


e 2 Hn{iyx)] =dn{l + 0{-) 


1 


c 


n / V — 1 


1/2 




where 




In order to use this formula for calculating Pn,tun, we use Eq. ([3]) changing 
the integration variable to get 


/ oo / s 2 

= 4a„7ri.e-"'/2j,(3.^-i)/3(i + 0(1))4 


where 


In = 


i: 


x^ — 1 


1/2 


Ai2(z2^/3c)dx. 


An easy power expansion leads to 


4a„7ri.e-"V2j,(3.^-i)/3 = „i/3 + 0(1)^ , 


therefore 


Pn,tun = 2"/V/3 (^1 + 0(1)^ 4. 

It is convenient now to introduce a new integration variable 

t = v^l\. 

Eig. [3] shows the dependence of C on x. We hnd 


dC f x^ — 1 

dx V C 


therefore 


and thus 


1/2 


roo 

In = / fn{t)A[‘^{t)dt, 


(14) 
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Figure 3: Plot of the function ((x) for x > 1. 


where 


fn{t) 


C 

— 1 


is implicitly defined through C = v and Eq. (I12p that allows us to find x 
as a function of Since 


.-4/3 = (2n + l)-2/3 = 2-2/3n-2/3 (^1 + 0(i)^ , 


for Pn,tun we get 


P„.„n = 


(15) 

Denote by Fn the integral in Eq. (|15p: 


Fn = 

poo 

/ /„(t)Ai2(t)dt. 

Jo 

(16) 

We will show now that 



~ Foo = 2-2/3 

1 

Jo “ 22/332/3r(l)2’ 

(17) 
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and therefore 


Pr 


n,tun 


32/3r(i)2 


n 


-1/3 


in agreement with Eq. dS]). Fig. 0] shows numerically calculated ratios 



Figure 4: Numerically calculated ratios Foo/Fn, 6 < n < 500. 


Foo/Fn for 6 < n < 500. In order to verify Eq. (fT7)l . we notice that (see the 
Appendix) the function 


f{x) 


(|x(x2 - 1)V2 _ 3posh 

— -^, yx > 1 ) 


(18) 


is monotonically decreasing, from /(I) = lim 2 .^i_|_ f{x) = 2“2/3 0 at 

X — )• oo. It follows that the sequence fn{t) is uniformly bounded. Moreover, 
with fixed t > 0, and since C = (2n + when n —oo, C ^ 0, thus 

X —>■ 1, therefore lim„^oo/n(^) = /(I) = Thus 

( lim fn{t)] Ai^(t)dt. 

\n—)-cxD / 

by the dominated convergence theorem. 


lim 


Fn = 


f 
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4 Next order approximation 

Using the method described in the previous section one can easily get the 
next term in the asymptotic formula for F^. The result is 

/-oo n-2/3 -2/3 /-oo -i 

Fn = 2~‘^/'^ / Ai^{t)dt --- / tAi^{t)dt + 0{—), (19) 

Jo 5 Jo ™ 

which, after some simplihcations leads to the formula 

Pn,tun = ^ (0.133975 - + 0 {n-P^). ( 20 ) 

The second term in Eq. ()19ll can be derived from the second term in the 
Taylor series expansion /„(t) = /n(0) + {dfn/dt)\t=o t + if we use Eq. 

(HI and notice that 

{dfn/dt)\t=o = i^"'‘^^(d/n/dC)lc=o ~ 

Fig. [5] shows the comparison of the exact data with the two curves. The 
black curve, obtained from Eq. (l20]l is closer to the data than the grey 
one, obtained from the simple approximation given by Eq. ([8|). Additional 
correction terms have been obtained by Paris in m by different means. 

5 Appendix 

Fig. [6] shows the computer generated plot of the function Q{x)/{x‘^ — 1) for 
X > 1. While the behavior of this function can be easily guessed from the 
plot, it is not that easy to prove rigorously that the function has the prop¬ 
erties that are necessary for the application of the dominated convergence 
theorem used in this paper. 

Lemma 1. Let f{x) be as in Eq. m8\) . Then f is monotonically decreasing, 
from /(I) = lim 3 ;^i+ f{x) = 2~‘^P to 0 at x ^ oo. 

Proof. We have 

fix) = > (21) 

where 

~ x(x^ — — cosh“^(x) 

ii^) = (x2 - 1)3/2 ’ 
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n 

Figure 5; Details of the comparison for n > 512 : exact data (dots), simple 
formula -Pn,tun = 0.133975(upper curve), and the formula with the 
first correction term Pn,tun = 0.133975— 0.0122518(lower curve). 



Figure 6: Plot of the function C,{x)/{x^ — 1). 
and it is enough to show that the statement in the lemma holds for /. Using 
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the identity 


cosh ^(x) = log(x + \/x2 — 1), 

and introducing s, 0 < s < 7r/2, through x = secs, we can write /(x) as a 
function z(s) of s 


z(s) 


sec s tan s — log(sec s + tan s) 
tan^ s 


a(s) 

b(s)’ 


( 22 ) 


and to prove the lemma we need to prove that z is monotonically decreasing 
from z(0) = lims^o+ z(s). Let 


h(u) = secutan^ u. 


By differentiation of the formulas below it is easy to verify that 

h(u)du, b(s) = 3 / h(u) secudu. 

Jo 

We notice that a(0) = b{0), and that, by de PHospital’s rule 



z(0) 


a'(s) 


i,// \ 

s—>-0+ b (s) 


2 

3 sec 0 


2/3. 


Therefore, using Eqs. (fT8)) . (|2T]l . /(I) = 2 Taking the derivative of Eq. 
(f2^ we get 

_ a'{s)b{s) - a{s)b'{s) 

b^s) 

The denominator is positive for s > 0. The numerator is negative owing to 
the fact that, for s > 0, 


a'{s)b{s) — a{s)b'(s) = 6 ^h{s) J h{u) sec u du — h{s) sec s j h{u)di 
= dh{s) / /i(u)(sec u — sec s)du, 

Jo 

with h{s) > 0, and sec(u) < secs for u < s. Therefore z\s) < 0. 
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